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We describe a new approach to computing energy lev- 
els of a non-relativistic bound-state of two constituents with 
masses M and m, by a systematic expansion in powers of 
m/M. After discussing the method, we demonstrate its po- 
tential with an example of the radiative recoil corrections to 
the Lamb shift and hyperfine splitting relevant for the hydro- 
gen, muonic hydrogen, and muonium. A discrepancy between 
two previous calculations of 0(a(Za) 5 m 2 /M) radiative recoil 
corrections to the Lamb shift is resolved and several new terms 
of 0{a{Za) 5 m 4 /M 3 ) and higher are obtained. 



The theory of non-relativistic bound states in QED 
remains an important source of information about fun- 
damental physical parameters, like the fine structure con- 
stant and the masses of the electron, muon and proton, 
among many others Simple atoms, which are being 
studied in laboratories, differ significantly in the ratios of 
their constituent masses. Two situations can be distin- 
guished. The first one is the case when the masses of the 
two constituents of the bound state are equal, with the 
positronium as the most important example. The sec- 
ond case is a bound state with two very different masses, 
e.g. hydrogen, muonium, muonic hydrogen. Both sit- 
uations represent two special limits of a general mass 
ratio case. In both limits certain simplifications are pos- 
sible. In the context of this Letter, the case of equal 
constituent masses was discussed to some extent in ^J,|) . 
In this Letter we consider the case when the masses of 
the constituents differ significantly from one another. 

Our main goal is a practical algorithm which allows 
a calculation of the bound-state energy levels in a given 
order of perturbation theory (in a and Za) as an expan- 
sion in powers and logarithms of m/M with an arbitrary 
precision. The opposite situation, i.e. calculation of the 
energy levels to all orders in a but in a fixed order in the 
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ratio m/M, has been studied in the literature 

In many practical situations only the first few terms 
of the expansion in m/M are required. Nevertheless, we 
believe that it is useful to construct such an algorithm in 
its generality. First, higher corrections in the ratio m/M 
might become relevant. For example, in the muonic hy- 
drogen m/M corresponds to m ll /M v ~ 0.113, not a very 
small parameter. In exotic hadronic atoms, such as pionic 
hydrogen, that ratio might be even larger. In hydrogen 
and muonium, where m/M is smaller, the very high pre- 
cision of experiments warrants a precise computation of 
the recoil effects. Second, a complete algorithm means 
that one can obtain the whole series in m/M at once and 
additional cross checks become possible. 

We first recall that non-relativistic bound-state ener- 
gies can be computed using an effective field theory ||. 
In we have shown how dimensional regularization 
facilitates this approach. One has to distinguish two dif- 
ferent contributions. The first one is the contribution of 
the relativistic region in the loop-momentum integrals; in 
what follows we will refer to these contributions as "hard" 
contributions. This is usually obtained as a Taylor ex- 
pansion of the relevant scattering amplitudes in spatial 
momentum components of external particles (which can 
be taken on-shell). 

Second, there is the so-called "soft" contribution, given 
by usual time-independent perturbation theory in quan- 
tum mechanics. An important point to note is that the 
soft contribution can in general be easily evaluated for 
arbitrary masses of constituents. As one can see from 
the Schrodinger equation, this is so because the essential 
dynamics of a non-relativistic bound state is governed 
by the reduced mass of the system rather than by the 
masses of individual constituents. In this respect, for the 
soft contributions the relation between the two masses 
is not very important and once the equal mass case has 
been solved, the rest follows easily. 

Therefore, in the situation where the two masses are 
different, the real problem is in computing the hard con- 
tribution and this is what we are going to discuss in this 
Letter. We will show that there is a simple way to ex- 
pand the hard scattering diagrams in powers of m/M. 
The essential advantage of this method is that it can be 



automated and many terms of the expansion can be easily 
computed. The only limitation is the available computer 
power; as a matter of principle, infinitely many terms in 
the m/M expansion can be obtained. High-performance 
symbolic algebra software is of great help in such com- 
putations (we use FORM @). 

The remainder of this Letter is organized as follows. 
First the method is described in detail. Next, we com- 
pute the a(Za) 5 radiative recoil corrections to both the 
Lamb shift and the hyperfine splitting up to the fourth 
order in the expansion in m/M. Finally, we present our 
conclusions. 

The method we are going to discuss is based entirely 
on using dimensional regularization. Note that for consis- 
tency one also needs the soft contribution in dimensional 
regularization; as we mentioned earlier this part of the 
problem is well understood [||,[|. 

The hard diagrams should be evaluated exactly at the 
threshold (zero relative velocity of the constituents). For 
this reason, the relevant Feynman integrals depend on 
only two scales, m and M. Since we are interested in 
their expansion in m/M , it is useful to be able to expand 
the integrands, so that only homogeneous one-scale in- 
tegrals have to be evaluated. Once this is achieved, the 
calculations simplify dramatically. 
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FIG. 1. The forward-scattering radiative-recoil diagrams. 
The bold line represents the heavy constituent of the 
bound-state (e.g. proton if we consider hydrogen) and the 
thin line — the light one (an electron). Diagrams with the 
crossed photons in the f-channel are not displayed. 

Our method is motivated by the known procedure 
which permits an expansion of Feynman diagrams in 
large momenta and masses 

Although originally 
formulated in a different way, that procedure can be re- 
formulated more practically using the notion of momen- 
tum regions. To arrive at the result one should follow 
a sequence of steps (1) determine large and small 

scales in the problem; (2) divide the entire integration 
volume into regions where each loop momentum is of the 
order of some of the characteristic scales; (3) in every re- 
gion perform a Taylor expansion in the parameters which 



are small in the given region; (4) after the expansion, ig- 
nore all the constraints on the regions and perform the 
integration over the entire integration volume; (5) add 
the contributions of different regions to obtain the final 
result. The only step in this sequence which might ap- 
pear counter-intuitive is the step 4, since one may sus- 
pect some double counting. The reason why that does 
not happen is that the scale-less integrals vanish in the 
dimensional regularization. This in turn implies that the 
results obtained from the integrals over different regions 
are different analytic functions of the parameters of the 
problem. Below this procedure will be demonstrated in 
some detail. 

To illustrate the method we focus on the last diagram 
in Fig. |l| and consider the following scalar integral: 
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Here [d D k] stands for d D k/(2ir) D , Pl = mQ, p 2 = MQ, 
where Q = (1,0,0,0) is the time-like unit vector. Only 
the relevant infinitesimal imaginary parts of the propaga- 
tors have been displayed. We are going to illustrate the 
expansion of the integral in Eq. (Q) in powers of m/M 
following the five steps outlined above. 

There are four momentum regions to be considered. 
In the first one all the momenta are of the order of the 
large mass M . In this case one can expand the electron 
propagators in mQki. The resulting integrals are all of 
the form: 
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with some integer powers a^. One immediately recognizes 
that all these integrals are identical with the general two- 
loop self-energy integrals of the particle with the mass M 
for which the general solution is known Jlifl . 

Next, there are two momentum regions where either 
ki ~ M and k 2 ~ m or k\ ~ M and k 2 ~ M, but 
k\ + k 2 ~ to. It is then easy to see that after a Taylor 
expansion in the small variables, the integral factorizes 
into a product of two simple one-loop integrals. 

The fourth region is determined by the condition k\ ~ 
k 2 ~ m. In this case the heavy particle propagator can 
be expanded in powers of k\ and in essence it becomes a 
static propagator. The general integral in this case has 
the form 
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Such integrals represent the only new type required for 
this calculation and the easiest way to solve them is to 



employ the integration-by-parts techniques p.5| , [t6| . Any 
integral J can be algebraically expressed as a combina- 
tion of the two-loop on-shell self-energy integrals and four 
new master integrals. The latter are the only integrals 
we have to compute, but this can be easily accomplished 
with help of Fcynman parameters. The results read 
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Let us add that all momentum regions other than the 
ones we discussed lead to scale-less integrals and are 
therefore not relevant. This concludes the construction 
of our expansion algorithm. 

We have applied this algorithm to compute the 
0(a(Za) 5 ) radiative recoil corrections to the Lamb shift 
and the hyperfme splitting of a general QED bound state 
composed of two spin- 1/2 particles with the masses m 
and M. It is well known that in this case the soft contri- 
bution is absent and the hard corrections shown in Fig. 1 
are the only diagrams we have to consider. We have done 
the calculation in a general covariant gauge; the cancella- 
tion of the gauge parameter dependence serves as a check 
of the computation. 

For the 5-wave ground state energy E we define: 
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where J = 0, 1 is the total spin of the two fermions form- 
ing the bound state. 

For the hyperfme splitting we obtain: 
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where /z = mAf /(m + M) is the reduced mass of the 
bound state. 

For the spin-independent energy shift we find 
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To our knowledge the terms C(to 3 /M 3 ) and higher 
are new for both E^fs and E aver , while the other terms 
have been obtained previously. In addition, the coeffi- 
cient of the 0(m/M) term in Eq. (j^) has been subject 
of some controversy, since two different numerical results 
have been reported, |l7|-|l9|] and |2(J . 

Our result for this term, 
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is in excellent agreement with the numerical result of 
Ref. HjJ where the coefficient -1.324029(2) was ob- 
tained. 



The discrepancy in the Q(m 2 /Ma(Za) 5 
to the Lamb shift reported in 



corrections 
19 1 and EG] has been 



the major source of the theoretical uncertainty in the 
so-called isotope shift (apart from the uncertainty asso- 
ciated with the proton and deuteron charge radii, see 
below), i.e. the difference between energies of 25 to 15 
transitions in deuterium and hydrogen: 

AE = [£(25) - E(1S)] D - [£(25) - E(1S)] H . (10) 



Experimentally, AE is known with the uncertainty of 
about 0.15 kHz pi[| ; the theoretical uncertainty associ- 
ated with higher order QED effects and with the un- 
certainties in the electron-to-proton and the electron-to- 
deuteron mass ratios is about 1 kHz each. On the other- 
hand, the difference in the results of Refs. ]l7|-|l9[ and 
pof leads to 2.7 kHz difference in AE. Our result for 
this term, Eq. removes this discrepancy in favor of 
the result of Ref. §§. 

It is well known that the high accuracy of the exper- 
imental value of AE cannot be used directly because of 
significant uncertainties in the value of the proton and 
deuteron charge radii which enter the theoretical for- 
mula for AE. In this situation the problem is usually 
turned around and one determines the difference of the 
charge radii of the proton and deutron using AE. Here, 
we do not pursue this topic any further. The related 
phenomenology can be extracted from Ref. [pifl ( sec also 
the recent review [ p2| , where results of Ref. |20(| should 
be used). 

We have constructed an efficient algorithm which per- 
mits an expansion of the energy levels of a bound state 
of two constituents with masses m and M in powers of 
m/M. This expansion is similar to, although not iden- 
tical with, the asymptotic expansions of Feynman dia- 
grams familiar from particle physics. We have demon- 
strated the usefulness of this procedure by computing 
several terms in the m/M expansion for the a(Za) 5 ra- 
diative recoil corrections to both the Lamb shift and the 
hyperfine splitting of a general QED bound state. 

Although we have only described a calculation of the 
radiative recoil corrections, the method is clearly appli- 
cable to all other types of corrections relevant for the 
bound states. In particular, the pure recoil corrections 
can be treated in a similar way. It remains to work out 
the details in that case, but the principles are clear. 

One of the terms in our result for the radiative recoil 
corrections to the Lamb shift is the 0(a(Za) 5 /j, 3 /(mM)) 
term for which two different numerical results have been 
previously reported. Our calculation confirms the result 
of Ref. @. 

Another aspect of this work might be related to higher 
number of loops. It is clear that the described method 
can be systematically applied in higher orders. Probably 
more important, it may facilitate the extraction of terms 
enhanced by In M/m which can be determined from the 
singularities of the contributions of different expansion 
regimes. Since those singularities must cancel in the com- 
plete result, their coefficients can be found by a partial 
calculation of the divergent parts of those contributions 
which can be evaluated most easily. 
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